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Abstract 

Stabilizing unstable periodic orbits in a chaotic invariant set not only reveals information about its struc- 
ture but also leads to various interesting applications. For the successful application of a chaos control 
scheme, convergence speed is of crucial importance. Here we present a predictive feedback chaos control 
method that adapts a control parameter online to yield optimal asymptotic convergence speed. We study 
the adaptive control map both analytically and numerically and prove that it converges at least linearly to a 
value determined by the spectral radius of the control map at the periodic orbit to be stabilized. The method 
is easy to implement algorithmically and may find applications for adaptive online control of biological and 
engineering systems. 

PACS numbers: 05.45.Gg, 02.30.Yy 



*Electronic address: bick@nld.ds.mpg.de 

I 



1 



I. INTRODUCTION 



Some chaotic attractors contain infinitely many unstable periodic orbits. These can be seen as a 
"skeleton" for the chaotic attractor, therefore revealing important information about the dynamics 
of the system itself. By suitable perturbations the stability of these unstable periodic points can 
be changed; a control perturbation renders them stable. Such "chaos control" has applications in 
many fields [1], including biological [2] and artificial neural networks EH). 

In the last twenty years, different methods for stabilizing unstable periodic orbits have been 
suggested. The seminal work by Ott, Grebogi, and Yorke (OGY) [5] and its implementations 
employ arbitrary small perturbations of a parameter of the system to stabilize a known unstable 
periodic orbit of a discrete time dynamical system. A successful application of the OGY method, 
however, requires prior knowledge about or online analysis of the dynamics to determine fixed 
points and their stability properties. 

A different approach is given by predictive feedback control (PFC) (6113 which overcomes this 
disadvantage. In this approach the future state of the dynamics calculated from the current state is 
fed back into the system to stabilize a periodic orbit. This feedback control is noninvasive, i.e., the 
control strength vanishes upon convergence, and is extremely easy to implement. It is a special 
case of a recent effort to stabilize all periodic points of a discrete time dynamical system J8l HI 
which is also closely related to nonlinear successive over-relaxation methods ifTOlfTTTl . It has been 
extensively studied Ifl2l - rr6ll and extended [fT71[T8l with respect to its original purpose as a tool for 
examining the structure of chaotic attractors. 

In any real world application, speed of convergence is of crucial importance. For example, if 
a robot is controlled by stabilizing periodic orbits in a chaotic attractor [3], the time it needs to 
react to a changing environment is bounded by the time the system needs to converge to a periodic 
point of a given period. Hence, in praxis, one desires to tune the control parameter such that the 
spectral radius of the unknown periodic point which the system converges to is minimized. To the 
best of our knowledge, previous works on chaos control have not considered convergence speed 
while maintaining its simplicity in terms of implementation. Adaptation of the control parameter 
has an impact on convergence speed. However, existing adaptation mechanisms [|3l [191 have two 
major shortcomings; they do not optimize for speed and, for adaptation of heuristic nature, may 
adapt the parameter to regimes where stabilization fails. 

Here we introduce an adaptation method that overcomes these shortcomings. It adaptively 
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tunes the control parameter online to achieve optimal asymptotic convergence speed. This work is 
organized as follows. In the second section we review the PFC method and introduce the notation 
that will be used throughout the paper. In the third section, we present the adaptation method and 
prove its convergence properties. As an example, the well-known logistic map is studied both 
analytically and numerically in Sections 4 and 5 before giving some concluding remarks. 

II. PRELIMINARIES 

A differentiable map / : W 1 — > IR n gives rise to a dynamical system through the evolution 
equation 

Xk+l = f(Xk) (1) 

with %t E W 1 for all k E Z. The sequence (xf.), k E N is called an orbit of the dynamical system 
with initial condition x and if f op (xk) = Xk+ P for all k > we say that the orbit is periodic with 
period p. Here, 

r = /o/o-o/ 

S s, ' 

p times 

denotes the p-fold composition of /. Let Fix(/) := {x E R n | f(x) = x} be the set of fixed 
points, i.e., periodic points of period one. Note that any periodic orbit is a fixed point of the p-th 
iterate of the map / so we will use the expressions fixed point and periodic orbit interchangeably. 

Let A C IR 71 be a forward invariant subset of M n with respect to /, i.e., f(A) C A. If periodic 
points are dense in A and / maps transitively, then we call A a chaotic set. Julia sets [|20| , as 
described below, are examples of such chaotic sets. Let df\ x denote the derivative of / at x E M. n 
and id : IR n — > W n the identity map. 

The results of [9] are now summarized as follows: 

Proposition II.l. Suppose x* E Fix(/) and the derivative df\ x * and df\ x * —id are real, nonsingu- 
lar and diagonalizable. Then there exist parameters fx > and an orthogonal matrix E 0(n) 
such that x* is an attractive fixed point of the map obtained from f through the transformation 

S(/J,,C k ) : / ^ id+p,C k (f - id) =: g^. 

In particular, S(p, Ck) preserves the set of fixed points, that is Fix(/) = Fix(<7 M ). 

In fact, it can be shown that the number of matrices Ck needed to stabilize all fixed points of 
a given map / is quite limited. They depend on the local stability properties of fixed points and 
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Control parameter p 



Figure 1: Sketch of the dependence of the spectral radius p x * (p) on p for some fixed point x* according to 



Proposition II. 1 



there are types of fixed points that can be stabilized for Ck G {± id}. For a given x* G Fix(/) 
with Xj E C, j G {1, . . . , n}, being the eigenvalues of dg^\ x * we want to denote by 

p x 4fi):= max {\Xj\} 

J6{l,...,n} 

the spectral radius, i.e., the maximum of the absolute values of the eigenvalues of the derivative of 
at x*. We have 

dg fl \ x = id +/iC fc (d/J x — id). (2) 

for all x G W 1 . In other words, the proposition above ensures the existence of p, and Ck(x*) 
for a given x* G Fix(/) such that the transformation S(p>,Ck) gives p x *{p) < 1, cf. Figure [Tj 
Therefore, with these parameters, the fixed point x* of / is an attracting fixed point for g M . 

The results above are directly related to predictive feedback chaos control methods. A transfor- 
mation T v : / i— > g is called a chaos control transformation if g can be written as g = f + r\c with 
control perturbation c : E n — > M. n and rj G K. Note that in case Ck G {± id} the transformations 
S(C k ,ii) are chaos control transformations since 

M = / + (l=F^)(id-/) 

with 77 = 1 — /i. Therefore, we will refer to these transformations S(p:, Ck) as PFC transformations. 
Without loss of generality we restrict ourselves to the case Ck = id. 

The results of Proposition |II. l[ however, give little information about the speed of convergence, 
except for the fact that when decreasing /1 towards zero convergence takes longer and longer as 
the spectral radius approaches one. In the vicinity of a stabilized fixed point, convergence is at 



least linear and the rate of convergence is bounded from above by the quantity p x * (/i). In order to 
obtain an adaptation method that increases the speed of convergence we therefore have to minimize 
p x * (p) using the control parameter p. For a random initial condition, we do not know to what fixed 
point x* (if any) the trajectory will converge to. We only have a converging sequence Xk — > x*. In 
other words, we are looking for a way to obtain a sequence pk — > Poo where 



Hoa = SUp < p > 



and assumptions of 

Vfi>0:p x .(fi)<p x *(fi) } (3) 



Prop. II. 1 are satisfied 



the optimal p to minimize p x *(p). Define Aoo := p x *(p 00 ). 

In applications, the control parameter p plays a double role; on the one hand, it can be used to 
turn chaos control on and off, p = 1, on the other hand it is the crucial parameter to stabilize the 
periodic orbits and to determine the speed of convergence. 

Define the class of functions 

F(p ,p) := {/ | card ({x G Fix(f p ) I pM < 1 for C fc = id}) > 0} , 

with parameters p G N, and p > and where card denotes the cardinality of a set. The sets 
^(pojp) are the functions / with a chaotic set that have at least one periodic orbit of period p 
which can be stabilized for the given parameters. 



III. AN ADAPTATION METHOD TO ACCELERATE CHAOS CONTROL 

In this section, suppose / G J(/io,p) for some p > and without loss of generality, p = 1, 
since we can replace / with the p-th iterate. Suppose is the transformed map after applying 
S(p, id). Furthermore we assume that for all times k < the system evolves according to Equa- 
tion ([T]), i.e., with r] = 1 — p = 0, along a trajectory of points in the chaotic set A. At time k = 
the control parameter p is set to po. Therefore, because of the assumptions on /, there is at least 
one periodic orbit of period p on the chaotic attractor which is now an attracting periodic orbit. 
Let sFix(/) denote the set of these stabilized fixed points. 
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A. Close to a fixed point 



Recall two facts: any differentiable map h E on an open set [/ C 1" is Lipschitz- 

continuous on any compact K C U, that is 

\\h(x) - h(y)\\ < \\dh\\ K -\\x-y\\, 

for all x, y E K where || • ||^ is the supremum of the operator norms induced by a norm || • || on 
K and dh is the total derivative. Furthermore, for any contraction h on a Banach space (X, || ■ ||), 
i.e., a map that satisfies 

\\h(x) -h(y)\\ <L\\x-y\\ 

with a Lipschitz constant L < 1, the Banach Fixed Point Theorem gives the existence of a unique 
fixed point x* together with the error estimates — Xk\\ < ^tz \\ x o ~ x i\\ an d \\xk+i — Xk\\ < 
L \\xk — Xk-i\\. Here, Xk = h ok (xo) for an initial condition xq E X. 



Let x* E sFix(/) be fixed. According to Proposition II. 1 there exists a Ao < 1 for /i = fi 
sufficiently small such that p x *(fio) < \q. Therefore, there exists a vector norm || • || such that 
we have HdgJ^* || op < A for the induced operator norm cf. Figure [2| a). We will omit the index 
indicating the operator norm when it is clear from the context. Henceforth all norms denote this 
vector norm and the induced operator norm respectively. 



Let B(e, x) denote a ball of radius s centered at x and by B(e, x) its closure. Assume that for 
e small enough there is a constant K > such that 

ll|d^|,.||-||d^| x .|||<K||x*-x|| (4) 



for all x with \\x* — x\\ < e independent of \i. This condition is depicted in Figure |2[b). Now we 
can choose e < e such that Hdg^ UH gT^i) < 1- Put differently, for sufficiently small So > there 
exists an e E (0, e) such that Wdg^ {xWw^- ?y < A + 5o =: L < 1. The choice of e (corresponding 
to the size of the ball around x*) depends on A , /U , and S . 

Remark III.l. Condition (|4]) is satisfied in case f has a bounded second derivative on sFix(/) 
due to the functional dependence of the derivative of on /i as given by ([2]). 

Algorithm III.2. For given f E J-^/icp), /io, A , and K let Xq E B(e,x*). The convergence 
acceleration algorithm consists of the following steps: 

Step 1 (Iterate): Calculate X\ = g w (x ). 




Figure 2: (a) An appropriately chosen norm approximates the spectral radius from above with p x *{pLo) < 
||d<7/xola:*|| — Ao- (b) According to Condition Q, Hdg^l^ || lies in a K \\x* — Xfc||-tube around Hdf?^* ||. 
While iterating, this tube becomes smaller and smaller. As another consequence of Q, we have 
Hdg^yi g^ in a Ke-tube around Hd^l^* ||. Therefore, for e small enough, there is Lq < 1 such that 

ll d 0W)U II B(e,z-) - L °" 

Step 2 (Optimize fi): Minimize the " cost function" ||d^| ;Cl || with respect to fi G (0, 1) under the 
conditions 

l{p) : = ||d# M U| + (gg) ||x - xi|| < L (5a) 
p, maximal (5b) 

where L = A + #o above, cf. Figure \3\ 

Step 3 (Set quantities): If the minimization under constraints of Step 2 returned a result fi opt 
then set ^ := /i opt ,X 1 : = Wdg^ || + ||s -^i||. S x := (jz^) INo-^ill and 

Li := Ai + 5\. Otherwise set jii := /i , Ai := A , S x := S and Li := L . 

Repeat the steps with all indices increased by one. 

For this method we obtain the following results. 



Lemma III.3. With the assumptions of Algorithm III. 2 above we have that for any initial condition 



xq G B(e, x*) Algorithm III.2 yields a trajectory Xk — > x*. 



7 




Figure 3: (a) Inequality ( pa) does not always have to be satisfied, (b) Since — > x* we have that after a 
finite time N the function l{fj) is below Lq for some {i. 

Proof. If the optimization process does not give a result, convergence is ensured by Proposi- 
tion jTLT and the Banach Fixed Point Theorem. Without loss of generality, suppose that optimiza- 



tion yields a result for k = 1. Then because of Q and (|5a]) we have 



\\^9^\x4 bE^%^) < ll d £«U*ll +K\\x 1 -x*\\ 
< Hdfl^JaJI + 2K - x*\ 



( 2KL \ 
\l-L J 



< Hd^Jxill + z z~ - xo 



= Li < L < 1. 

Therefore, x 1 is contained in a ball around the fixed point x* on which the map g Ml is a contraction 
with contraction coefficient L\. The same calculation is valid for subsequent optimization steps 
for k > 1. □ 

The lemma above ensures that the adaptation does not compromise convergence against the 
stabilized fixed point. But will optimization actually take place? For a map with K = 0, adaptation 
is not necessary since || dg^ \ x * \\ = || dg^ \ Xk || and therefore we can set \i straight to the optimal value. 



Lemma III.4. With the assumptions of Algorithm \III.2\ Inequality ( |5a| ) z'^ satisfied every finitely 
many steps. 

Proof. By definition we have Hd^l^* || < L = \ + So with So > 0. Hence we have 
11^5^0 U* II < ^o- Let ( > be such that £ < L — Hd^l^. ||. Since \\x k — x k+1 \\ is a Cauchy 



sequence and Hdg^JaJI — > \\dg^ \ x *\\ there is an N G N such that \\%n — a^iv— 1 1| < § and 



lll d 0/J*'ll - lld^ol^HI < §. Thus, we have 



Therefore, Inequality ( |5a| ) will be satisfied after maximally N = N steps. 

Inductively, by increasing all indices above by N, the same argument gives a sequence iVj, 
/ G N, of indices for which Inequality ( |5al ) is satisfied. This completes the proof of the assertion. 

□ 

Remark III.5. Although the adaptation method gives a sequence p,f- that minimizes the norm 
while ensuring convergence, it is not clear how often optimization yields a result. Additional 
conditions on the map f, such as requiring monotonicity of Hdg^Uj. || in Xk, influence how often 
the parameter p, will be adapted. On the other hand, additional constraints make the theory less 
broadly applicable. 



If Inequality ( paj ) is satisfied for some k > 0, then, because of continuity, it holds for a whole 
closed neighborhood of /i k . This gives p, k+1 with Hd^^J^H < Hd^J^H unless Hd^l^H is 
constant on that interval. 

Definition III.6 (ED). A matrix norm \\ ■ \\ on R nxn is called minimal for M £ R nxn if p(M) = 
\\M\\. 

The main results of this section can now be summarized in the following theorem. 

Theorem IIL7. Suppose f G J-"(/io,p) for p > such that f satisfies Q with K > in a 
neighborhood of x* G sFix(/). Furthermore, lete,Xo, and 5q be chosen as described above. Then 



for any initial condition xo G B(e, x*) Algorithm 111.2 minimizes an upper bound for the spectral 
radius p x * (/i). 

In particular, if the induced operator norm \\ ■ \\ op is minimal for dg^ \ x *, it converges at least 
linearly with asymptotic convergence speed A^. 

Remark III.8. For dimension n = 1, the Euclidean norm is minimal. 



Proof, (of Theorem III.7 ) Lemmas III.3 and III.4 ensure convergence against the fixed point x* 



and adaptation of the control parameter ji after a maximum of some finite number of steps. 
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By construction, jik tends to a value which minimizes the norm of the derivative of g M at x*. 
For arbitrary dimension n we have p x *(pk) < H^SWls*!!- ^ m addition the norm is minimal, 
in the limit the spectral radius is minimized yielding optimal asymptotic convergence speed, i.e., 



Remark IIL9. One could also use convergence acceleration transformations R22\l in order to get 
a better approximation to p x *(p,). However, to exploit the acceleration within the framework of 
this theory, one would have to have suitable error estimates for the transformed sequence. 

The choice of the size neighborhood B(e, x*) depends on the desired estimate of the contraction 
constant. It is clearly bounded from above since we have to make sure that there is a contraction. 
On the other hand, it is desirable to take a neighborhood as large as possible to make the method 
applicable to as many initial conditions as possible. 

B. From local to global 

We want to consider the situation where the control is turned on at a random point in time. We 
choose indices such that this time is k = 0. In general, the initial condition xq for the adaptation 
method is unknown and it is likely to be outside of a neighborhood B(e,x*) as defined in the 
section above. One possible scenario is the existence of a chaotic attractor that makes up part 
of the old chaotic attractor and ending up in its basin of attraction when the control parameter is 
turned up. Another likely scenario is to have x close to the boundary of the basin of attraction of 
one of the stabilized fixed points. An initial condition close to the basin boundary implies a long 
transient iteration before the adaptation method becomes applicable. 

We want to quantify the latter scenario. Since we assumed the system to follow the evolution 
equation Xk = f{xk~i) for all k < 0, xq is distributed on the attractor according to some f- 
invariant measure m on A. Suppose m is an ergodic probability measure on A. So m(U) is the 
probability that xq e U at time k = for any measurable subset U C A. 

Let B(e(x*), x*) denote the neighborhoods of x* 6 sFix(/) for which the acceleration method 
described in the previous section is applicable for all initial conditions within that neighborhood. 
Because of f £ ^(/io, p) at least one of these balls is not empty. Define 



Pk -> At. 



'OO- 



□ 



V 
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to be the part of the union of all these neighborhoods on the attractor. Thus, if V is measurable, 
m(V) is a lower bound for the probability that the adaptation method described above converges 
if the parameter fx is set to fx at a random point in time and dynamics evolved on the chaotic set 
before. Furthermore, we define Vk := \Ji<k9^o (V). Now Pk = m(Vk) is a lower bound on the 
probability that the algorithm will converge after letting the transformed system evolve for k time 
steps after being initialized with /i = /i at time k = 0. 

As k — > oo tends to infinity the set Vk will converge to the union of the basins of attraction of 
the stabilized fixed points. Hence, we obtain a function 

ip(fi ) := lim m(14) 

depending on the initial parameter /i . The value lim inf^^o vKa'o) = f° r some (p G [0, 1] 
determines the size of the basin of attraction of the stabilized fixed point. 



IV. ADAPTIVE PFC FOR THE LOGISTIC MAP 



As an example, we apply the PFC transformation to the logistic map given by the quadratic 
polynomial £ r (x) = rx(l — x) with the real parameter 1 < r < 4. It is well known that there 
are parameter values for which the dynamics is chaotic on some subset A of the unit interval 
/ = [0, 1]. In particular, for r = 4 the whole unit interval is a chaotic set. Here, we study the 
period one orbits; higher periods can be treated similarly. 



A. Calculating the Adaptation Parameters 



First, we want to calculate the quantities for the adaptation method described in ^ III A The 
second derivative of t r exists everywhere on IR and is bounded on compact subsets. Within this 
section let h! denote the derivative of a differentiable function h. Here, we treat the cases Ck G 
{± id} simultaneously by allowing the control parameter fx for the transformed function 

g^, r {x) = S(ii,C k )(£ r )(x) = x + ix{l r {x) - x) 

to range within the interval [—1,1]. For \i = 1 we obtain the original system and around fi = 
either, where Ck = id when ji is positive and Ck = — id when \i is negative. 
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Since \g" (x)\ = |^' r '(x)| = 2r \n\ for all x £ I, the maximum in p, is taken for = 1. 
Therefore, if we set K = 8, we obtain a constant independent also of the parameter r and the sign 
ofC fc . 

The two fixed points of / r are x* = and x* = The derivatives at the fixed points are 
g'^ r (0) = 1 + p(r — 1) and g'^i^) = 1 — p(r — 1). Hence, x* = is stable for p negative 
(Ck = — id) and x* = — for /i positive (C^ = id). To apply the adaptive method, the initial 
parameters need to be determined as in §111 A for a given p the bound Aq can be calculated 



directly from the derivative. Furthermore, we have to find e that defines a neighborhood of x* 
for the initial condition x and a given initial p . From the local stability and Q we obtain that 
convergence is ensured if 

A = 1 + |/i | (r- 1) > -1, (6a) 
ATe- |/i | (r- 1) < 0, (6b) 
Ke+|/i |(r-l) <2 (6c) 

For either x* this gives \fj,o\ < This results in a bound for the size of the neighborhood of x* 
in which the map is a contraction, 

. f 2-/i (r-l) 
s<mrn| , — ^- 

The optimal bound £ < is achieved for /i (r — 1) = 1. 

It is desirable to choose e as large as possible (to cover as many initial conditions as possible) 



while keeping the whole expression on the left hand side of ( |6a| ) inequality as small as possible (a 
smaller contraction constant L leads to stronger contraction). This choice depends on the initial 
guess (J®. 

Remark IV. 1. The chaotic set A depends on the choice of the parameter r, so we obtain a family 
of chaotic sets A r . Note that we do not necessarily have G A r or E A r . We have A± = I 
so that the two fixed points are contained in A*. Otherwise, for a given fixed point x*, e has to be 
chosen large enough such that A r C\ B(e,x*) ^ 0. 

The constant parameters K — 8 and /x , Ao, and e as given by ([6]) together with an approxima- 
tion of the measure on A are the basis of the calculations of lower bounds for the probability of 
convergence in the following section. 
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B. PFC on the complex plane 



The quadratic polynomial defining the logistic map can also be seen as a polynomial over the 
complex numbers. Iteration of complex polynomials is a classical example in complex analytic 
dynamics and the theory developed there can tell us something about the effect of the PFC trans- 
formation S(fx, Ck) for Ck G {± id}. Here, this geometric point of view allows us to calculate the 
full basin of attraction of the stabilized fixed points. In particular, we also obtain convergence for 
general, complex valued initial conditions in a neighborhood of the periodic orbits in the complex 
plane. 

Recall some notions from one-dimensional complex dynamics l|2~0ll . Suppose / : C — > C is 
holomorphic. A point z G C = CU{oo}is said to be in the Fatou set F(f) if there is an open 
neighborhood of z on which the family of iterates {f ok \ k G N} is normal. Its complement is 
called the Julia set J(f) and constitutes the boundary of all Fatou components that contain any 
stable periodic orbits. Both of these sets are forward and backward invariant with respect to the 
map /. The Julia set is a chaotic set in our sense. Henceforth, we denote the complex variable by 
z. 

Let / G C[z) be a complex polynomial. Note that the result of the PFC transformation g^ r again 
is a polynomial of the same degree in the complex variable z unless / is constant or ji = 0. The 
logistic map is defined as a polynomial of degree two. The dynamics of quadratic polynomials 
are conjugate to the dynamics of a polynomial z 2 + c where c is a complex parameter. The 
parameter c can be characterized in terms of the orbit of the only finite critical point z = (since 
(z 2 + c)'(0) = 0), where the points for which that orbit is bounded constitute the Mandelbrot set 
M. For c G M. there can be bounded Fatou components corresponding to the basin of attraction 
of a stable periodic orbit. 

The logistic family described above is conjugate to the subset [—2, |1 of the intersection of M. 
with the real axis. Since g^ r again is a real quadratic polynomial of degree two for £ r and these 
polynomials only keep the real axis invariant if c 6 E, every g^ r is conjugated to a quadratic 
polynomial z 2 + c with real c. For given r, the relationship between this complex parameter c and 
the control parameter /i is given by 

c r (/i) = i(l-/i 2 (r-l) 2 ) 

for /i / 0. Hence, varying the parameter fi results in a "shift" of the dynamics up the real axis 
until it approaches c = | as fi — > 0. From the equation above, one can also see that the dynamics 
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Figure 4: Julia sets for parameters fi = —0.65 (left) and /i = —0.2 (right) for Ck = — id. The Julia set 
°f 67t,4 is depicted in dark gray whereas the Julia set of the original map, i.e., the unit interval, is depicted 
in medium gray. The fixed points z* = and z* = | are marked by black dots. Shaded circles indicate 

B(e, z*) for [1q = fi. 

of are conjugated for fj, — 1 and fj, — — 1, the former case corresponding to the unperturbed 
system. 

What does stabilization of fixed points mean in terms of complex analytic dynamics? An 
unstable fixed point is contained in the Julia set. The goal of stabilization is to turn this fixed 
point into a stable one, i.e., that it now belongs to a bounded Fatou component. That is, the 
transformation should deform the Julia set in a way that it does not contain the targeted periodic 
point anymore. 

Let us consider the case r = 4 in more detail. The Julia set J(£ r=4 ) = / is equal to the whole 
unit interval. The probability distribution m is given by a beta distribution with parameters both 
equal to | (cf. for example Il23l0 . i.e., with probability density function 

p(x) = ^7TX2(1 — x)2 j 

Suppose Ck = id and fi small enough. In this case z* = | is the stabilized fixed point. In the 
previous section we calculated the maximum size of the ball around the fixed point for which the 



14 





Re(z) 



1 





Re(z) 



1 



Figure 5: Julia sets for parameters fi = —0.65 (left) and fx = —0.2 (right) for Ck = id. The Julia set of 
g Mj 4 is depicted in dark gray whereas the Julia set of the original map, i.e., the unit interval, is depicted 
in medium gray. The fixed points z* = and z* = | are marked by black dots. Shaded circles indicate 

B(e, z*) for //q = fj,. 

adaptation method works straight away. This radius is given by e < | and 



P = m (V ) < m(B(e, z*) f) I) = m 



5 7' 



dx 



| 7ra;2 (1 — x) 2 



7T 



arctan 



vf) ~ arccot (yvj 



0.1895. 



In Figure HI one can see that the whole unit interval is contained in the bounded Fatou component. 
Backward iteration takes this set closer to the boundary of this Fatou component. This means that, 
for /i small, 

<pM = 1, 

Pq < Pf. — m(Vjt) < 1 and <p — 1. Therefore, trajectories will converge to the stabilized periodic 
point with probability one. 

The picture is slightly different for Ck = — id and ji small enough. Now, z* = is the stabilized 
fixed point. Again we have e < | and therefore 



P = m (y ) < m{B{e, z*) n /) = m 
2 



0. 



dx 



7T£2 (1 — x) 2 



= -arccot ~ 0.2301. 

In this case backward iteration yields a different result as can be seen in Figure |5j Part of the set 
of initial conditions A is in the basin of attraction of infinity and the intersection with the Julia 
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set is exactly the fixed point z* = f . Therefore, the probability to converge with a random initial 
condition is less than one. Integrating the probability density function gives 



3 

4 ax 2 



'0 7T£2 (1 — x) 2 3 

Therefore, we have (p = | and P < m(Vk) < |. In contrast to the case Ck = id this means that 
for Ck = id and z* = a trajectory with an initial condition on / distributed according to m will 
diverge with a probability of one third. 

When considering higher periods of such a polynomial map, the Julia sets are more complicated 
as the degree of the iterated polynomial rises exponentially. The situation changes qualitatively 
when considering the predictive feedback control dynamics of higher dimensional maps by inter- 
preting them as functions C n — > C n . In general, the dynamics of holomorphic, higher-dimensional 
maps is more diverse since even low-dimensional invertible maps give rise to complicated dynam- 
ics El. 



V. NUMERICAL RESULTS 

To compare the speed of the adaptive method (ACC) with the original PFC chaos control in a 
real world application, we performed numerical simulations for the logistic map £4. The results 
for Ck = + id, /i € [0, 1] and periods one and two are summarized in Figure |6j One can clearly 
see that for most initial values of the control parameter, the adaptive method yields an increase 
in convergence speed. The results for Ck = — id and period one are similar but the convergence 
probability is lower (not shown) in accordance with the results of the previous section, cf. Figure[5} 
In the case of period p = 2, the orbits stabilized by negative ji are the period one orbits, which is 
reflected in our numerical results (not shown). A non-optimized, ad hoc choice of parameters for 
the adaptive method of K = 8 and L Q = 0.99 (independent of the initial condition) were employed 
in the simulations. The criterion for convergence time T was given by \xt — xt-i\ < 1CT 10 but 
reliability was determined after checking for the correct period. 

The convergence reliability, i.e., the fraction of trials where the above criterion is fulfilled after 
some time T, is not improved by the adaptive method. However, it is possible to amend the adap- 
tation method to lead to convergence for most initial conditions xq within the convergent regime, 
independent of the initial value of the control parameter ^0 (the modified method is denoted by 



ACCD). When adapting, the ACC method has to check whether Criterion pa) is fulfilled. If this 
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Figure 6: Speed and reliability comparison of the original PFC chaos control (// = fixed) and both 
ACC and ACCD for the logistic map £4 with p = 1 (left) and p = 2 (right). Times are only plotted 
if more than 1% of initial values that lead to convergence to the correct period. Gray shading indicate 
the convergence reliability; dark gray corresponds to all methods converging, light gray to the reliability of 
ACCD. Convergence time T is given by \xt — %t-i\ < 10 -10 calculated for 1000 random initial conditions 
after a transient of random length. 

is not the case after M iterations, the modified method simply reduces fi to a certain fraction v. 
To prevent \i of becoming too small, we imply a threshold 9 below which /i cannot decay. Put in 
other words, the modified method ACCD will automatically decrease ji towards zero to reach the 



convergence regime if Inequality ( |5a| ) is not satisfied within a given number of steps. 

The modified method ACCD behaves like the original ACC method for initial values of fx in 
the convergent regime while leading to convergence outside of it, cf. Figure [6] (here M = 50, 
v = 0.7, 9 = 0.1 for period p = 1, and 9 = 0.05 for period p = 2). Failure of convergence that 
is due to the existence of a range of diverging initial conditions, however, will persist, even with 
the decay. The results are similar for a broad parameter range (e.g., decay rate v E [0.65, 0.99] and 
decay kick-in time M E [10, 100]). For a decay rate too close to 100% or a too large decay kick-in 
time, it will take many iterations to reach the convergent interval. On the other hand, if the decay 
kick-in time is too small or does not exist at all, decreases even if it is in the convergent interval 



as Criterion (5a) is not fulfilled all the time, unnecessarily increasing convergence time. 
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VI. DISCUSSION 



Here we presented a method which adapts the control parameter p, of the PFC method in or- 
der to accelerate the systems convergence to a periodic orbit. In contrast to heuristically chosen 
methods, our adaptation does not compromise convergence. The algorithm converges in a neigh- 
borhood of every periodic orbit that was stabilized by the stabilizing transformation. Assuming 
the existence of an invariant, ergodic probability measure on the chaotic attractor, we obtain an 
analytic bound for the probability P k that the system converges to a periodic orbit if the chaos 
control is switched on at an arbitrary point in time. Although these results are stated in the frame- 
work of discrete time dynamical systems, they can also be applied to stabilize continuous time 
systems after discretization such as taking Poincare sections. The logistic map provides an exam- 
ple for which we can calculate the parameters for the method. We estimated the probability of 
convergence and highlighted its dependence on the fixed point to be stabilized and the associated 
matrix C^. 

Our method was stated in the general context of "chaotic sets." In general, such sets do not need 
to be local or even global attractors of the dynamical system. In fact, the Julia sets considered 
in the example are repelling rather than attracting. In applications, however, an attractor would 
be desirable such that the process of stabilization becomes repeatable. That is, after the control 
perturbation is turned off by choosing the appropriate value for the control parameter, the dynamics 
would return to the attractor and the process can be started over again. 



Apart from its role as a chaos control method, the estimates described in \ III B give information 
on the PFC method itself. It allowed us to calculate the size of the basin of attraction for varying fi 
in our example. Decreasing ji always leads to slower convergence since the eigenvalues converge 
to one as p, — > 0. So is it possible to find an optimal p for a given map? Since any adaptation 
method increases the computational cost of the chaos control method, a priori estimates of such 
crucial quantities are of importance. Furthermore, the choice of the stabilization matrix Ck de- 
pends on the type of fixed points in the chaotic attractor. Hence, global statistics for a given map 
/ of the periodic orbits and their stability properties might yield some a priori estimates. 

Our numerical studies suggest that it is possible to get reliable convergence without a priori 
knowledge of the exact values for the parameters. A slight modification of the method yields a 
hybrid method that finds the regime of control parameter in which the dynamics converge online 
and then adapt it to the optimal value. This simplification, however, comes at a cost in conver- 
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gence speed. By definition, PFC cannot distinguish between a periodic orbit of period p and any 
q\p, a divisor of p. Our numerical calculations, however, indicate that this does not influence re- 
liability of the chaos control method. This is most likely caused by the exponential growth of the 
number of periodic orbits. In the future, it would be desirable to add a mechanism that rigorously 
distinguishes between the target period and its divisors to prove optimal convergence. 

An adaptation method for chaos control is a step towards solving the intuitively contradictory 
problem of optimizing speed while maintaining simplicity in implementation. However, as dis- 
cussed above, it leads to further challenging research questions that have to be addressed in the 
future. 
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